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Abstract 

An important class of problems exhibits smooth behaviour in space 
and time on a macroscopic scale, while only a microscopic evolution law 
is known. For such time-dependent multi-scale problems, an "equation- 
free framework" has been proposed, of which the gap-tooth scheme 
is an essential component. The gap-tooth scheme is designed to ap- 
proximate a time-stepper for an unavailable macroscopic equation in 
a macroscopic domain; it uses appropriately initialized simulations of 
the available microscopic model in a number of small boxes, which 
cover only a fraction of the domain. We analyze the convergence of 
this scheme for a parabolic homogenization problem with non-linear 
reaction. In this case, the microscopic model is a partial differen- 
tial equation with rapidly oscillating coefficients, while the unknown 
macroscopic model is approximated by the homogenized equation. We 
show that our method approximates a finite difference scheme of ar- 
bitrary (even) order for the homogenized equation when we appropri- 
ately constrain the microscopic problem in the boxes. We illustrate 
this theoretical result with numerical tests on several model problems. 
We also demonstrate that it is possible to obtain a convergent scheme 
without constraining the microscopic code, by introducing buffer re- 
gions around the computational boxes. 



1 Introduction 

For an important class of multi-scale problems, a separation of scales ex- 
ists between the (microscopic, detailed) level of description of the available 
model, and the (macroscopic, continuum) level at which one would like to 
observe the system. Consider, for example, a kinetic Monte Carlo model 
of bacterial growth [23] • A stochastic model describes the probability of an 
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individual bacterium to run or "tumble" , based on the rotation of its flag- 
ellae. Technically, it would be possible to evolve the detailed model for all 
space and time, and observe the macroscopic variables of interest, but this 
would be prohibitively expensive. It is known, however, that, under certain 
conditions, one can write a closed deterministic equation for the evolution of 
the macroscopic observable (here bacteria concentration, the zeroth moment 
of the evolving distribution) as a function of macroscopic space and time. 

The recently proposed equation-free framework JIJ can then be used 
instead of stochastic time integration in the entire space-time domain. This 
framework is built around the central idea of a coarse time-stepper, which is 
a time-Ai map from coarse variables to coarse variables. It consists of the 
following steps: (1) lifting, i.e. the creation of appropriate initial conditions 
for the microscopic model; (2) evolution, using the microscopic model and 
(possibly) some constraints; and (3) restriction, i.e. the projection of the 
detailed solution to the macroscopic "observation" variables. This coarse 
time-stepper can subsequently be used as "input" for a host time-stepper 
based algorithms performing macroscopic numerical analysis tasks. These 
incude, for example, time-stepper based bifurcation code to perform bifur- 
cation analysis for the unavailable macroscopic equation [23 EH El EEEj ■ A 
coarse timestepper can also be used in conjunction with a projective inte- 
gration method to increase efficiency of time- integration |Z|. This approach 
has already been used in several applications [2S1 E| > and also allows to do 
other system level tasks, such as control and optimzization |24j . 

When dealing with systems that would be described by (in our case, 
unavailable) partial differential equations, one can also reduce the spatial 
complexity. For systems with one space dimension, the gap-tooth scheme 
|14j was proposed; it can be directly generalized in several space dimensions. 
A number of small intervals, separated by large gaps, are introduced; they 
qualitatively correspond to mesh points for a traditional, continuum solution 
of the unavailable equation. In higher space dimensions, these intervals 
would become boxes around the coarse mesh points, a term that we will also 
use throughout this paper. We construct a coarse time-Ai map as follows. 
We first choose a number of macroscopic grid points. Then, we choose a 
small interval around each grid point; initialize the fine scale, microscopic 
solver within each interval consistently with the macroscopic initial condition 
profiles; and provide each box with appropriate (as we will see, to some 
extent artificial) boundary conditions. Here, we constrain the macroscopic 
gradient to a value that is determined by the macroscopic solution field. 
Subsequently, we use the microscopic model in each interval to simulate 
until time At, and obtain macroscopic information (e.g. by computing the 
average density in each box) at time At. This amounts to a coarse time-Ai 
map; this procedure is then repeated. 

This "coarse" scheme has been used with lattice-Boltzmann simulations 
of the Fitzhugh-Nagumo dynamics |131 114j and with particle-based simula- 
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tions of the viscous Burgers equation |Hj . It was analyzed in the case where 
both the microscopic and the macroscopic model are pure diffusion |14j . 
where it was shown to be equivalent to a standard finite difference scheme 
of order 2 in space, combined with an explicit Euler step in time. Here, we 
extend the analysis for the gap-tooth scheme in several ways. We derive a 
formulation which approximates difference schemes that have higher order 
accuracy in space; and we analyze the convergence of this generalized scheme 
for a one-dimensional parabolic homogenization problem with non-linear re- 
action. In this case, the microscopic model is a partial differential equation 
with rapidly oscillating coefficients. The macroscopic model is the effective 
equation that describes the evolution of the average behaviour. In the limit, 
when the period of the oscillations becomes zero, this effective equation is 
the classical homogenized equation. The goal of the gap-tooth scheme is to 
approximate the effective equation by using only the microscopic problem 
inside the small boxes. We analyze the accuracy of the method analytically 
for the case where the homogenized solution is close to the effective solu- 
tion. This analysis is important, because it shows that the gap-tooth scheme 
approximates the correct effective equation in the presence of microscopic 
scales. 

It is worth mentioning that many numerical schemes have been devised 
for the homogenization problem. Hou and Wu developed the multi-scale fi- 
nite element method that uses special basis functions to capture the correct 
microscopic behaviour |10[ lll|. Schwab, Matache and Babuska have devised 
a generalized FEM method based on a two-scale finite element space |22U19| . 
Runborg et al. [20] proposed a time-stepper based method that obtains the 
effective behaviour through short bursts of detailed simulations appropri- 
ately averaged over many shifted initial conditions. The simulations were 
performed over the whole domain, but the notion of effective behaviour is 
identical. The guiding principle in equation-free timestepper-based compu- 
tation is to perform numerical tasks on an unavailable equation. The time 
derivative for the evolution of the field is not obtained from a formula, but 
estimated from observations of short, appropriately initialized and processed 
detailed dynamic simulations in (portions of) space. When more informa- 
tion about the structure of the unavailable equation is known (e.g. that it 
is a conservation law for a known observable), it makes sense to modify the 
general time-stepper based procedure appropriately; one can, for example, 
estimate the time derivative based on flux computations using an available 
microscopic simulator in (portions of) space. This modification of equation- 
free computations for the case of conservation laws forms the basis of the 
generalized Godunov scheme of E and Enguist jS] and of the finite difference 
heterogeneous multiscale method of Abdulle and E pQ. Our approach has 
focused on the general case where the structure of the unavailable equation is 
not known. It is interesting to pose the question about how one might know 
whether a system is effectively a conservation law (and additional ques- 
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tions, such as whether a system is Hamiltonian, or, possibly, integrable). 
A computer-assisted methodology for the equation-free exploration of such 
questions is introduced in [T5] . 

In the gap-tooth scheme discussed here, the microscopic computations 
are performed without assuming such a form for the "right-hand-side" of the 
unavailable macroscopic equation; we evolve the detailed model in a subset 
of the domain, and try to recover macroscopic information by interpolation 
in space and extrapolation in time. We note again that the gap-tooth scheme 
as it is presented here, is only a part of the equation- free solution framework. 
In this paper we examine the properties of this coarse time-stepper per se; 
yet one should keep in mind that the coarse time-Ai map will eventually be 
used inside a projective integration code, or a bifurcation/continuation code. 
The combination of gap-tooth timestepping with projective integration has 
been termed patch dynamics |14j . 

This paper is organized as follows. In section |2J we discuss a general 
order formulation of the gap-tooth scheme. Subsequently, in section |21 we 
discuss some basic theoretical results on mathematical homogenization, and 
we give a relation between the averaged solution and the homogenized so- 
lution. In section 0J we analyze the convergence of our method for the 
model homogenization problem. Numerical results confirming the theorem 
are shown in section El This section also contains some examples for which 
the theory is strictly speaking not valid. We discuss a modified version of 
the gap-tooth scheme in section H3 that avoids constraining the macroscopic 
gradient during simulation. We introduce so-called buffer regions that shield 
the dynamics inside each box from boundary effects. At the outer boundary 
of the buffer box, one can subsequently apply whatever boundary conditions 
the microscopic code allows. We propose to study the resulting scheme by 
the (equation- free) numerical computation of its damping factors. We show 
how this can be done for a diffusion problem with Dirichlet boundary con- 
ditions. We conclude in section where we also point out some next steps 
of this research. 



2 The gap-tooth scheme 

We consider a general reaction-convection-diffusion equation with a depen- 
dence on a small parameter e, 

d ( x \ 

—u t (x,t) = f lu e (x,t), —u e (x,t), -^u e (x,t),x,-j , (1) 

with initial condition u e (x,0) = uq(x) and Dirichlet boundary conditions 
u e (0,i) = v and u e (l,t) = v±. We further assume that / is 1-periodic in 

y = f- 
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We are only interested in the macroscopic (averaged) behavior u(x,t), 
which is a "filtered" version of u e (x,t). To this end, we define an averaging 
operator for u e (x,t) as follows, 

U(x, t) := S h (u)(x, t) = T ' « e (£. t)d£. (2) 

This operator replaces the unknown function by its local average in a small 
box of size h >> e around each point. If h is sufficiently small, U(x,t) 
should be a reasonable approximation to u(x,t). 

The averaged solution U(x, t) satisfies an (unknown) evolution law, which 
we assume also diffusive, 

^U(x, t)=F (u(x, t), ^-U(x, t), ^U(x, t),x; h^j . (3) 

Note that this equation depends on the box width h. 

The goal of the gap-tooth scheme is to approximate the solution U (x, t), 
while only making use of the detailed model Q. Suppose we want to obtain 
the solution of the unknown equation (jHJ) on the interval [0,1], using an 
equidistant macroscopic mesh II(Ax) : = {0 = Xo < X\ = xq + Ax < . . . < 
xn = 1}- We construct a time At-map for U(x,t) in the following way. 
We consider a small box (tooth) of length h « Ax centered around each 
mesh point, and solve the original problem in each box. To determine 
the simulation within each box completely, we impose boundary constraints 
and an initial condition as follows. 



Boundary constraints. Each box should provide information on the evo- 
lution of the global problem at that location in space. It is therefore crucial 
that the (artificially imposed) boundary conditions are chosen to emulate the 
correct behaviour in a larger domain. Since the microscopic model (^Q) is dif- 
fusive, it makes sense (thinking of traditional explicit numerical schemes) to 
impose a fixed macroscopic concentration gradient at the boundary of each 
small box during a time interval of length At. We determine the value of 
this gradient by an approximation of the macroscopic concentration profile 
u(x, t) by a polynomial, based on the (given) box averages i = 1, . . . , N. 

u(x,t n ) zipi(x,t n ), x e [Xi - ~,Xi + -], 

where (x, t n ) denotes a polynomial of (even) degree k. We require that the 
approximating polynomial has the same box averages as the initial condition 
in box i and in | boxes to the left and to the right. This gives us 

***** Pi(C, t n )dC = U? +j , j = - k -... k - (4) 

i+j 2 
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One can easily check that 



S h ( P 1){x,t n ) = £ U^L^x), L{,-(x) = f[ ,[ X _T\ ^ 
3=-| ! =-i 

where - {x) denotes a Lagrange polynomial of degree k. The derivative of 
this approximating polynomial is subsequently used to obtain the value of 
the gradient at the boundary of the box. 

If we did have an equation for the macroscopic behaviour, we would use 
these slopes as Neumann boundary conditions. Here, we use these deriva- 
tives to constrain the average gradient of the detailed solution u(x, t) in box 
i over one small-scale period around the end points, 

1 p^-f+f 8 1 P'+f+f 8 

s _ w u«,m = .,-, -J ^m-*. m 

x ' 2 2 X '~^ 2 2 

Note that we approximate a box average in a box of macroscopic size 
h >> e, while we average for boundary condition purposes over a length 
scale e that is characteristic for the microscopic model. Hence, we replace 
each boundary condition and its effect on the simulation by an algebraic 
constraint. 



Initial conditions. For the time integration, we must impose an initial 
condition u l {x, t n ) in each box [xj — 4, Xj+§], at time t n . We require u % (x, t n ) 
to satisfy the boundary conditions and the given box average. We choose a 
quadratic polynomial u(x,t n ), centered around the coarse mesh point Xj, 

u l (x,t n ) = a(x - Xi) 2 + b(x - Xi) + c. (8) 

Using the constraints (JJJ) in the limit for e — > and requiring 

1 



h 



It — h. 



we obtain 



a= s A 8 ±. h= s i +s i c = Tj n - —(s + - s~) (9) 
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Figure 1: A schematic representation of a gap- tooth time-step. We choose a 
number of boxes of size h around each macroscopic mesh point Xi and interpolate 
the initial averages (dots) in a number of boxes around Xi (dashed profile). The 
derivatives at the boundary and the average are used to create an initial profile in 
box i (full line). 



The algorithm. The complete gap-tooth algorithm to proceed from U n 
to U n+1 is given below (see also figure 0: 

1. Lifting. At time t n , construct the initial condition u*(x,i n ), i = 
0, . . . , N, using the box averages U™ (j = 0, . . . , N) as defined in @. 

2. Evolution. Compute u l (x,t) by solving the equation Q until time 
t n+ i = t + At with the boundary constraints 0. 

3. Restriction. Compute the box average U™ +1 = j- f Xl+ h w e (£, t n+ i)d£ 
at time t n +±. 

It is clear that this amounts to a "coarse-to-coarse" time-At map. We write 
this map as follows, 

U n+1 = S k (U n ;t n + At), (10) 

where S represents the numerical time-stepping scheme for the macroscopic 
(coarse) variables and k denotes the degree of interpolation. 

Microscopic simulators. It is possible that the microscopic model is not 
a partial differential equation, but some microscopic simulator, e.g. kinetic 
Monte Carlo or molecular dynamics code. In fact, this is the case where 
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our method would be most useful. In this case, several complications arise. 
First of all, the choice of the box width h becomes important, since there will 
generally exist a trade-off between statistical accuracy (e.g. enough sampled 
particles) and spatial resolution. 

Second, the lifting step, i.e. the construction of box initial conditions, 
also becomes more involved. In general, the microscopic model will have 
many more degrees of freedom, the higher order moments of the evolving 
distribution. These will quickly become slaved to the governing moments 
(the ones where the lifting is conditioned upon), see e.g. ^1E!> but it might 
be better to do a constrained run before initialization to create mature initial 
conditions P^llS], 

Third, as already mentioned, imposing macroscopically inspired bound- 
ary conditions is non-trivial [TJI|- Moreover a given microscopic code may 
come with one of several "standard" microscopic boundary conditions. We 
will therefore examine the effect of incorporating simulations with such stan- 
dard boundary conditions in a gap-tooth context, provided we extend the 
simulation in a buffer region surrounding the computational "tooth" . The 
solution in the buffer is not used in the restriction step. This variant is 
examined more closely in section |HI 

Finally, even determining which and how many macroscopically inspired 
boundary conditions are needed, is a delicate issue. This is related with 
the order of the partial differential equation, i.e. the order of the highest 
spatial derivative. A systematic way to estimate this, without having the 
macroscopic equation, is given in |15j . 

3 Model homogenization problem 

Here, we review some basic results from homogenization theory. We note 
that we are interested in finding the effective behaviour of the solution. In 
our setup, we know that for sufficiently small e the effective behaviour is close 
to the homogenized behaviour, which is the limit of the solution for e — > 0. 
Since in some cases, the homogenized equation can be found analytically, 
we use this equation as our reference for the effective behaviour. 

3.1 Standard homogenization theory 

As a model problem, we consider the following parabolic partial differential 
equation, 

d d / / x \ d \ 

— u e (x, t) = — U^-j g^u e (x, t)\+ g(u e (x, t)), (11) 

with initial condition u e (x,0) = u°(x) and suitable boundary conditions. In 
this equation, a(y) = a (^) is periodic in y and e is a small parameter. 
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Consider equation 1)11(1 with Dirichlet boundary conditions u e (Q, t) = vq 
and u t (l,t) = V\. According to classical homogenization theory [3], the 
solution to ()11() can be written as an asymptotic expansion in e, 



oo 



u t (x, t) = u (x, t)+^e i (ui(x, -,t)j, (12) 

i=l 

where the functions Ui(x,y,t) = Ui(x,^-,t), i = 1,2,... are periodic in y. 
Here, uo(x,t) is the solution of the homogenized equation 

d d ( d \ 

—u {x,t) = — la*—u (x,t) J +g(u (x,t)) (13) 

with initial condition uq(x,0) = u°(x) and Dirichlet boundary conditions 
uo(0, t) = vq and uo{l,t) = v±; a* is the constant effective coefficient, given 
by 

a* = J a(y) (l - ^x(l/)) dy, (14) 
and x(y) is the periodic solution of 

the so-called cell problem. The solution of 1)15)1 is only defined up to an 
additive constant, so we impose the extra condition 

^1 

X (y)dy = 0. (16) 



From this cell problem, we can derive u±(x,y,t) = ^g^xiv)- 

These asymptotic expansions have been rigorously justified in the clas- 
sical book Under appropriate smoothness assumptions, one can obtain 
pointwise convergence of no to u t as e — > 0. Therefore, we can write 

\\u e (x,t) - u (x,t)\\ < C e, (17) 

where = ll/C^)!!^ = max^ |/(x)| denotes the oo-norm of /. Through- 

out this text, whenever we use ||-||, we mean the oo-norm. 

It is important to note that the gradient of u(x, t) is given by 

Odd 

— u e (x, t) = -g^ u o{x, t) + q-ux(x, y, t) + 0(e), (18) 

from which it is clear that the micro-scale fluctuations have a strong effect 
on the local detailed gradient. Nevertheless, since Ui(x,y,t) are periodic in 
y, the gradient of the homogenized solution can be approximated by the 
averaged gradient over one period e of the medium. The error is bounded 
by 



d 1 f x+ 2 d 

—uo(x,t)-- —u e (£,t)d(, 

OX € Jx-^ OX 



< C' e. (19) 
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3.2 Homogenization and averaging 

The gap-tooth scheme introduces an approximation on two levels. The 
scheme computes an approximation to the evolution of the averaged macro- 
scopic quantities instead of an approximation to the solution of the true 
homogenized solution. Before considering how well the gap-tooth scheme 
approximates this averaged behaviour, it might be of interest to show how 
the averaged behaviour approximates the homogenized solution. 

Lemma 3.1. Foru(x,t) sufficiently smooth, the averaged function 



U{x,t)=S h {u)(x,t) := - 



x 2 



u(£,t)d£ 



can be asymptotically expanded in h as follows, 

fh\ 21 

U{x,t) = u( 



, , ^fh\ M 1 d 21 



t) 



Em 



We omit the proof, but this can easily be checked using Maple. 

Using this lemma, we consider the homogenization problem of section 
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dt 



dx 



x\ d 



ol^fo *) = aZ \ a 7 -ZI u e( x > *) + 9(u e (x, t)) 



dx 



(20) 



In this case, we can bound the difference between the averaged solution 
U(x,t) and the homogenized solution uo(x,t) in the following way. 

Lemma 3.2. The difference between the homogenized solution uo(x,t) and 

x+ — 

the averaged solution U(x, t) = J h 2 u(£, t)d£ is bounded by 



\\U(x,t)-u (x,t)\\ <dh 2 + C 2 e. 

Proof. We first make use of the asymptotic expansion ([12)1 for u(x,t) and 
the triangle inequality, and subsequently of lemma 13.11 



\\U(x,t) - u (x,t) 



< 



X 2 



£+4 



u(£,t)d£ - u (x,t) 



h 
2 



' 2 



uo(£,t)d£ - u (x,t) 



+ 0(e 2 



< — 



24 
h 2 



dx" 2 



+ C 2 e 



< — 



max 

24 xe[o,i] 



ff 2 
dx 



■^u (x,t) 



+ C 2 e 
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This concludes the proof. □ 

This shows that the averaged solution is a good approximation of the 
homogenized solution for sufficiently small box width h. 

4 Convergence results 

To analyze the convergence of the gap-tooth scheme, we solve the detailed 
problem approximately in each box. Because h >> e, we can resort to 
the homogenized solution, and bound the error using equation (|17|) . It is 
important to note that we only use the homogenized equation for analysis 
purposes. We never make use of the homogenized equation in the imple- 
mentation. 

We first relate the gap-tooth time-stepper as constructed in section |2 
with a gap-tooth time-stepper for which the box problem is the homogenized 
equation with Neumann boundary conditions. 

Lemma 4.1. Consider the model equation, 

d d ( ( x\ d \ 

—u e (x,t) = — ( a —u t (x,t)j +g(u e (x,t)) , (21) 

where a(y) = a ( | ) is periodic in y and e is a small parameter, with initial 
condition u e (x,0) = u°(x) and boundary constraints 

1 /-^-l+t d 1 d 

- h e Q-MZ,t)dt = s-, - h e -u^m = si. (22) 

For e — > 0, this problem converges to the homogenized problem 
d d ( d \ 

—u (x,t) = — la*—u (x,t) j +g(u (x,t)) (23) 

with initial condition uq(x,0) = u°(x) and Neumann boundary conditions 

d 



dx u (x,t) 



h = sf, (24) 

X — X i "2" 



and the solution of (j21j) - (|22j) converges pointwise to the solution of (|23(1 - 
((23}, with the following error estimate 

||« e (x,t)-«o(x,t)|| <C 3 e. (25) 

This lemma can be checked using the two-scale convergence method [2] 
or formally by making use of asymptotic expansions [3]. 

Using this lemma, we now estimate the difference between a gap-tooth 
time-step using (|21H22|) and a gap-tooth time-step using the homogenized 
box problem (t2lTOl . 
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Lemma 4.2. Define U n+l = Sk(U n ,t n + At) as one gap-tooth time-step on 
the full problem 1)21(1 with the box constraints 1)22)1 . and U n+l = Sk(U n ,t n + 
At) as one gap-tooth time-step on the homogenized problem 1)23)1 with bound- 
ary conditions 1)24)1 . When U n = U n , we have 



n+l 



n+l 



m 



Proof. Denote the solutions of l)2l j) -(j22 ]) and ([23 )1 - (|24 |) . with the initial condi- 
tion Ui(x,t n ) determined by the lifting step, ©-© as u\(x,t) resp. u l Q {x,t). 
We can then write 



U 



n+l 



u: 



n+l 



«o(^'*n+l)dC 



'L i 



< 



max 



(u e (C, Wl) - &o(£, *n+l)) d £ 
l \u e (£,t n+1 ) - «o(£, Wl)|| 
< C 3 e 

Here, we bounded the average over the interval [xi — ^,X{ + \\ by the max- 
imum, and subsequently used lemma 14.11 This is valid since we assumed 
U n = U n . Therefore, the initial condition for both box problems is the 
same. 

This proves the lemma. □ 
The averaged solution U(x,t) satisfies a reaction-diffusion-like equation 



g(u(£,t))d$. 



(26) 



We denote a forward Euler /spatial finite difference approximation for (J2l 

as 

U n+1 = S k (U n ,t n + At), 

with k the order of accuracy of the spatial finite differences. The following 
theorem compares a gap-tooth time-step U n+1 = Sk(U n ,t n + At) with a 
finite difference time-step. 

Lemma 4.3. We denote a finite difference approximation of order k for the 
evolution ofU(x,t) as 

U n+1 = S k (U n ,t n + At), 

and one gap-tooth time-step with homogenized box problems \2°J\2J$ as U n+1 = 
Sk(U n ,t n + At). The exact solution of the homogenized equation is denoted 
byu(x,t). WhenU n = U n = S h (u){x,t n ), we have the following estimate 



U 



n+l 



u 



n+l 



< C 4 At 2 
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Proof. Denote the solution (j23[) -(|24 |) with the initial condition Ui(x,t n ) = 
a(x — Xi) 2 + b(x — Xi) + c, determined by the lifting step ©-©, as Uq(x, t). 

• We write an expansion for the solution of (|23|) - (|24[) using the method 
of separation of variables. The solution can be decomposed as follows 

u l Q (x,t) = u l (x,t n ) + v\x,t), 

where v l (x,t) is the solution of 

d ■ d ( d ■ \ 

of(x,t) = — U*—v l (x,t)j +a* .2a + g(v^{x,t)), 

with homogeneous initial condition and homogeneous Neumann bound- 
ary conditions. The last term is due to the spatial derivatives of the 
initial profile. We write v l (x, t) as a Fourier series with time-dependent 
coefficients, 

v l (x, t) = ^-i + cos ( k n( x ~%i)) + Yl Kit) sin (l n (x -Xi)), 

n=l n=0 

with k n = 2 ^ and l n = ( 2ra + 1 ) 7r _ The Fourier modes satisfy the 
homogeneous Neumann boundary conditions, and form a set of spatial 
basis functions for the solution. The time-dependent coefficients are 
given by 

9 r Xi+ 2 /9-rrn \ 

= T I V^U)cO S (^- Xi ))^ 

K(t) = \^^t)^(^^^-xi^ 

Each coefficient can be found by solving an ordinary differential equa- 
tion that is obtained by taking the time derivative and replacing the 
time derivative of the solution by the right-hand side of the partial 
differential equation. 

• We use this analytical solution to obtain an explicit formula for one 
time-step of the gap-tooth scheme. We write 

u: +1 = I r\ 2 m,t n +At)dt 

r +2 (4(£,i n )+^(e,t n +Ai))d£ 



h 



2 



up + fi&L+M, (27) 
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where the last step is due to the definition of u l (x,t n ) and the zero 
averages of sine and cosine. Thus, we only need to consider the coef- 
ficient ag(i) in what follows. 

For ao(t), we get the following ordinary differential equation 



+4 a ■ 

2 



with initial condition a l (t n ) = 0. This yields 



d_ 

df 



* at 

2 







4a* • a + - 



h 

' 2 



+ a* • 2ah + 



Xi+ 2 



J-i 2 



s(4(£>t))d£ 



where we could discard the first term because of the boundary condi- 
tions. 

The resulting formula for a gap-tooth time-step is 



I pt=t n +At rxi+^ 



t=tn 



g(u^,t))d^dt (28) 



• We now wish to connect one time-step of the gap-tooth scheme with 
one finite difference time-step for the equation for the averaged solution 
U(x,t). We first notice that 



s- 



2h 

1 ( dp l k (x,t) 
2h 



dx 



dpUx,t) 



dx 



2 S h (p l k )(x,t) 



?=-- 

2 
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Therefore, the second derivative of the fe-th order polynomial that 
interpolates the box averages is equal to 2a. Due to symmetry, the 
first and second derivatives of this polynomial in Xi are the standard 
finite difference approximation of order k of the function U(x,t n ). This 
can easily be verified with Maple. This leads to the following formula 
for one finite difference time-step of the equation for the averages 



U 



ra+1 



U? + 2a* • aAt + GfAt, 



(29) 



where 



/-in 



Xl 2 



with u(x,t) the exact solution of (|13|). 
We now estimate 



G?At 



I rt n +At rXi-, 2 

h 



< 



+ 



3> i 



1 



g(u(£,t n ))At 
t n +At 



g(u^,t))d(dt 
t n +At 



g(u(Z,t))dt)dt 
<7 (4 (£,t)))did£ 



The first term can be written as 

i 

h 



t n +At 



At 



g(u(Z,tn))At- 

g(u(Z,t n )) - 

u=u{^,t n ) 



1 

At 



1 f x * + 2 1 dg 
h 



h 2 du 

2 



du 
~dt 



g(u(tt))dt)d( 
Atd£ 



t tr. 



< At 

< C'At 2 , 

where we used a Taylor expansion for g and the chain rule. The second 
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term is estimated as 



1 
h 

< 



t n +At 



(s(u(f,t))-0 («&(£,*))) dtd£ 

1 f X i + l ftn+ At 



< L 



t n +At 



L(«(£,i)-t$(£,i))dtd£ 
(u(£, t„) - «°(£, t n )) + (t - t n )Cdtd£ 



< C"At 2 , 



where we used Lipschitz continuity of g and the fact that (due to 
lifting) 



u({,,t n )d£, 



4&tn)d£. 



This proves the lemma. 



We now have the following result. 



□ 



Theorem 4.4 (Local error). Define U n+1 = Sk{U n ,t n + At) as the result 
of one gap-tooth time-step for ill]) , and U n+1 = Sk(U n , t n + At) as a finite 
difference time-step for \2b]) . When \Jf = IJf = Sh(u)(xi,t n ) (the exact 
solution at (xi,t n ) ), the difference is bounded by 



Pi +1 -U^\\ <C 3 e + C 9 At z 
Proof. This follows immediately by combining lemmas 14 . 1 1 and |4 . Ml 

Therefore, we obtain the following error bound. 

Theorem 4.5. IfU n+l = Sk(U n ,t n + At) is a stable finite difference scheme 
for H3), then U n+1 = S k (U n ,t n + At). Moreover, if U° = U° = S h {u°){x), 
the error with respect to the homogenized solution uo(x,t) of hl'J\) is bounded 
by 

\\U? - u (xi, t n )\\ < dh 2 + C 2 e + C s ^ + C 9 At + C 5 Ax k 

Proof. We start by splitting the error, based on the origin of the error con- 
tributions, 

||c/™ - U (Xi,t n )\\ 

< \\U? - U?\\ + ||t/J* - U(xi,t n )\\ + \\U( Xi ,t n ) + UQ (xi,t n )\\ 

< \\U? -U?\\ + C 5 Ax k + C 6 At + dh 2 + C 2 e. 



rn+ll 



□ 
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The last two terms follow from standard finite difference theory and from 
lemma l3~Tl The first term merits further investigation, 



u n -u n \\ = 


\\s k (u n - 


-\t n „ 1 + At)-S k (U n -\t n - 


-1 + At) 


< 


\\s k (u n - 


-\t n ^ 1 + At)-S k (U n - 1 ,t n ^ 


-i + At) 


+ 


\\s k (u n - 


-\t n ^ + At)-S k {U n -\t n . 


-i + At) 


< 


\\A k (U n - 


~ l - U 11 - 1 ) II +C 3 e + C 4 At 2 , 





where the last line is due to lemma 14.41 and A k is the error amplification 
matrix for the forward Euler /spatial finite difference scheme. Therefore, 
stability of the finite difference scheme is necessary and sufficient to bound 
the errors from previous steps. By induction, we obtain 

\\U n -U n \\ < Cn(C 3 e + C A At 2 ) 

< C^- t (C 3 e + C 4 At 2 ) 

We prove the theorem by combining all terms. □ 

The result clearly shows the interplay between the different approxima- 
tions; we have an error due to the (macroscopic) finite difference scheme, 
an error that arises because we consider box averages, and an extra error 
due to the setup of the box problems in each time-step. The introduction 
of averaged Neumann boundary conditions generates an error which is inde- 
pendent of the time-step. We therefore have to make a trade-off between the 
accuracy that is gained by taking shorter time-steps and the accuracy that 
is lost because of more frequent reinitializations. This will also be shown in 
the numerical experiments. Projective integration [7j can help in reducing 
this error, since then less re-initializations are needed. 



5 Numerical results 

We show convergence of the gap-tooth method for a diffusion problem with 
a rapidly oscillating diffusion coefficient (section 15. lj) , a reaction-diffusion 
system (section 15.2(1 . and a system with a rough non-periodic (random) 
diffusion coefficient (section I5.3JI . 

5.1 Periodic diffusion coefficient without reaction 

Consider the following model problem, 

^u e (x,t) = ^- (a(^)^u e (x,t)\ , o(^) = l.l + sin(27r^) (30) 

with e = 1 • 10" 3 , x G [0,1], initial conditions u e (x, 0) = 1 - 4(x - 0.5) 2 , 
and Dirichlet boundary conditions u e (0,t) = u e (l,t) = 0. To solve the 
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0.2 0.4 x 0.6 0.8 1 0.595 0.6 0.605 



Figure 2: Left: Solution of equation (|30|l with an implicit Euler scheme with Sx = 1- 
10 -5 and St = 5-10~ 7 (full line), and a gap-tooth solution (dashed) with Ax = 0.05, 
At = 2.5 • 10~ 4 and h = 0.01, with internal finite differences as for the full problem. 
The boxes indicate the domains where the computations are done in the gap-tooth 
scheme. Right: A zoom shows the presence of the microscopic fluctuations. 



microscopic problem, we use a standard finite difference discretization in 
space and an implicit Euler time-step, with parameters Sx = 1 • 10~ 5 and 
St = 5 ■ 10 -7 . The corresponding homogenized equation is given by 

d ( d \ 

— ( a*— u Q (x, t)j, a* « 0.45825686. (31) 

With respect to the theoretical setup of section two additional approx- 
imations are made during the computations: the time integration inside 
each box is not exact; and we have to use numerical quadrature formulas 
to obtain the box average at each restriction. The resolution of the inter- 
nal time-stepper is such that these effects are negligible with respect to the 
other sources of error that we wish to study. In our code, we used the trape- 
zoidal rule as quadrature formula. Figure [2] shows the solution of (|30|) and 
the gap-tooth solution with Ax = 0.1, At = 1 • 10~ 3 and h = 0.01 at time 
t = 0.02. 



Difference with respect to finite differences. We first compare the 
results of the gap-tooth scheme for <|3U|) with those of a finite difference 
scheme for 1)3 1JI with the same coarse parameters. We use a spatial inter- 
polation/finite differences of order k = 2, with a coarse mesh of Ax = 0.1, 
resp. Ax = 0.05, for decreasing box sizes h = 0.04, 0.02, 0.01, 0.005, and for 
time-steps At = uAx 2 , with v = 0.1, 0.2, 0.4. The results are shown in table 
^for Ax = 0.1 and in table [21 for Ax = 0.05. They clearly show an OQi 2 ) 
decrease of the error initially, with a slow-down for smaller h, due to the 
additional 0{-^) term. We also see that the decrease of convergence speed 
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At = O.lAx 2 


At = 0.2Ax 2 


At = OAAx 2 




error 


ratio 


error 


ratio 


error 


ratio 


h = 0.04 


5.4189 • 10" 4 




5.3755 • 10~ 4 




5.3568 • 10~ 4 




h = 0.02 


1.4296 • 10" 4 


3.79 


1.3815 • 10~ 4 


3.89 


1.3584 • 10~ 4 


3.94 


h = 0.01 


4.3169 • 10~ 5 


3.32 


3.8297 • 10~ 5 


3.61 


3.5885 • 10~ 5 


3.79 


h = 0.005 


1.8221 • 10~ 5 


2.37 


1.3334 • 10~ 5 


2.87 


1.0896 • 10~ 5 


3.29 


Table 1: Difference between the gap-tooth scheme for ()30[) and a finite difference 
scheme for (|31|) for order k = 2 and Ax = 0.1 at time t = 2 • 10~ 2 . 




At = O.lAx 2 


At = 0.2Ax 2 


At = OAAx 2 




error 


ratio 


error 


ratio 


error 


ratio 


h = 0.04 


5.6378 • 10" 4 




5.5060 • 10" 4 




5.4275 • 10" 4 




h = 0.02 


1.7152 • 10" 4 


3.29 


1.5293 • 10" 4 


3.6 


1.5641 • 10~ 4 


3.79 


h = 0.01 


7.2618 • 10" 5 


2.36 


5.3027 • 10~ 5 


2.88 


4.3236 • 10" 5 


3.32 


h = 0.005 


4.7638 • 10~ 5 


1.52 


2.8043 • 10~ 5 


1.89 


1.8247 • 10~ 5 


2.37 



Table 2: Difference between the gap-tooth scheme for (|30|) and a finite difference 
scheme for l|31|) for order k = 2 and Ax = 0.05 at time t = 2 • 10 -2 . 



is affected by the time-step. The error decreases less rapidly for smaller 
is, due to the additional error in each restriction step. Also note a smaller 
decrease for Ax = 0.05, because in this case At is also smaller. Note that 
the difference with respect to finite differences does not depend dramatically 
on At for this example, due to the absence of a reaction term (see theorem 
14.4)1 . By comparing tables ^ and |21 we also see that the difference between 
the gap-tooth scheme and the finite difference scheme is independent of Ax, 
for fixed h and v. 

The 0(57 )-term. The evolution of the difference between the gap-tooth 
scheme and the finite difference scheme is shown in figure 03 We see that 
we start with a constant error at time t = due to the averaging of the 
initial condition. Note that this error is not important if one compares to 
the exact averaged solution U(x,t). It is an artifact of comparing to the 
homogenized equation instead of the effective equation. However, if uo(x,t) 
evolves according to (|3*Tj). U(x,t) = Sh(uo)(x,t) evolves according to the 
same equation. Therefore, we can eliminate the 0(h 2 ) term, by comparing 
to U(x,t) instead of uo(x,t). This allows us to show that the stagnation in 
tables ^ and [2] is really e-dependent. We compare the results of the gap- 
tooth scheme and the finite difference scheme of order 2 for Ax = 0.05 and 
h = 0.02 at time t = 2 • 10~ 2 . We first keep At = 1 • 10 -3 fixed, and vary 
e = l- 10~ 3 , 2 • 10~ 3 , 4 • 10~ 3 . Subsequently, we fix e = 1 • 10~ 3 , and vary 
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Figure 3: Difference between the gap-tooth scheme for (|30[) with k = 2, Ax = 0.05, 
At = 1 • 10~ 3 , h = 0.02, and a finite difference scheme with the same coarse 
parameters for time t = 0, 4 ■ 10" 3 , 8 ■ 10" 3 , 12 • 10~ 3 , 16 ■ 10" 3 , 2 ■ 10" 2 . 





error 


ratio 


e = 4 • 10" 3 
e = 2 • 10~ 3 
e = 1 • 10" 3 


3.0574 • 10" 4 
7.6239 • 10- 5 
1.9710 • 10" 5 


4.01 
3.87 





error 


ratio 


At = 0.5 • 10~ 3 
At = 1 • 10~ 3 
At = 2 • 10" 3 


3.9304 • 10" 5 
1.9710 • 10~ 5 
9.9133 • 10" 6 


1.99 
1.99 



Table 3: Left: Difference between the results of the gap-tooth scheme and the 
forward Euler/spatial finite difference scheme of order 2 for Ax = 0.05, At = 
1 • 10~ 3 , h = 0.02, at time t = 2 • 10~ 2 , for e = 1 • 10~ 3 , 2 • 10~ 3 , 4 • 10~ 3 after 
subtracting the 0(h 2 ) error in the initial data. Right: Difference for e = 1 • 10~ 3 , 
At varying. 



At = 0.5 • 10" 3 , 1 • 10" 3 , 2 • 10~ 3 . The results are shown in table US For this 
simple example, the error decreases quadratically with e, because the error 
constant of the 0(e) term is zero. If we combine the results from both tables, 

2 

we clearly see a decrease in error according to O(^) for this example. 

Error with respect to solution of the homogenized problem. We 

also show the error with respect to the exact solution of the homogenized 
problem. For this purpose, we compute the homogenized solution with a 
second order finite difference approximation in space with 5x = 1 • 10 -5 and 
implicit Euler time-steps with St = 5 ■ 10 -7 . The gap-tooth scheme is used 
with box width h = 0.005, Ax = 0.2, 0.1 and At = uAx k with v = 0.4 and 
order k = 2. We compare the gap-tooth and the finite difference solution to 
the exact solution for the homogenized equation at time t = 2 • 10~ 2 . The 
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(a) k = 2, Ax = 0.1 (b) k = 2, Ax = 0.05 

Figure 4: Error of gap-tooth (full line) with h = 0.005 and finite differences (dashed) 
with respect to an exact solution for the homogenized equation for Ax = 0.1 (left), 
Ax = 0.05 (right) and At = OAAx 2 . 



results are shown in figure 01 It is clear that the error is similar to that of 
the finite difference scheme. Note however that the errors will increase when 
the 0(h 2 ) and O(eAt) terms in the error become dominant. 

Higher order discretizations in space. We repeat the same experiment 
for a gap-tooth scheme of order k = 4, which we compare to a fourth order 
spatial finite difference approximation, with an explicit Euler time-step. As 
parameters, we choose Ax = 0.1, 0.05 and At = z/Ax 4 , with v = 0.4 
and h = 0.01. In order to view the 0(Ax 4 ) behaviour, we need to choose 
At correspondingly small, which will affect convergence through the 
term. The results are shown in figure 00 It is clear for Ax = 0.1 that the 
scheme approximates the fourth order scheme. In this case, the time-step 
At = 4 • 10" 5 . However, for Ax = 0.05, At = 2.5 • 10~ 6 , the error is already 
completely dominated by the O(^) term. 

5.2 Periodic diffusion with non-linear reaction 

As a second example, we consider the following reaction-diffusion equation 

§i UeM = ^ (^^M)) +u e (s,t)(i - ^^), (32) 

where a(|) = 1.1 + sin(27r|) as in section IB~T1 and b(x) = sin(27rx) + 1.2. 
This model can be interpreted as a one-species logistic growth model with 
macroscopically varying capacity b{x) and a rapidly oscillating diffusion co- 
efficient a(^). We choose periodic boundary conditions and a constant initial 
condition, u°(x) = 0.7. The corresponding homogenized problem is given 
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Figure 5: Error of the gap-tooth scheme (full line) and a finite difference scheme, 
both of order 4, with respect to the exact homogenized solution at times t = 4- 10 -3 , 
8 ■ 10~ 3 , 12 • lO" 3 , 16 • 10" 3 , 2 • 1CT 2 , where Ax = 0.1 (left), Ax = 0.05 (right), 
At = 0.4Ax 4 and h = 0.01. 




Figure 6: Reference solution (full line) and gap-tooth solution (dashed/circles) for 
the problem (EU at time t = 2 • 10~ 2 . 



by 

j^(z,t) = -J-- (a*|u e (*,i)) +« e (x,t)(l - ! ^^), (33) 

with a* w 0.45825686. Figure El shows the solution of (|32J), as well as the 
result of a gap-tooth simulation with parameters Ax = 0.05, At = 1 • 10~ 3 , 
h = 0.01 for t = 2 • 10 -2 . The reference solution was computed with second 
order spatial finite differences and an implicit Euler time-stepper, with Sx = 
1 • 10~ 5 and St = 5 ■ 10~ 7 . 

In contrast to the example in section 15.11 there will be also an O(At) 
difference between the gap-tooth solution and the finite difference scheme. 
To show this, we compare the gap-tooth scheme with Ax = 0.1 and At = 
uAx 2 with a finite difference scheme for the homogenized equation with the 
same coarse parameters for v = 0.05, 0.1, 0.2, 0.4. We did this for box width 
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/i = 0.005 


h = 0.01 




error 


ratio 


error 


ratio 


At = 4 • i(r 3 


1.3842 • 10" 4 




1.3929 • 10^ 4 




At = 2- i(r 3 


7.9135 • 10~ 5 


1.75 


7.9792 • 10~ 5 


1.75 


At = i • i(r 3 


5.1103 • 10~ 5 


1.55 


5.1496 • 10~ 5 


1.55 


At = 5 • 10" 4 


3.8014 • 10" 5 


1.34 


3.7959 • 10~ 5 


1.36 



Table 4: Difference between the gap-tooth scheme for (I30|) and a finite difference 
scheme for <|31[) for order k = 2 and Ax = 0.05 at time t = 2 • 10 -2 . 





h = 0.005 


h = 0.01 




error 


ratio 


error 


ratio 


At = 4 • 10~ 3 


1.0812 • 10" 4 




1.0982 • 10~ 4 




At = 2- 10-' 3 


4.9371 • 10" 5 


2.19 


5.1135 • 10~ 5 


2.15 


At = 1 • 10-' 3 


2.1459 • 10" 5 


2.3 


2.3411 • 10~ 5 


2.18 


At = 5 • 10~ 4 


8.1325 • 10" 5 


2.63 


1.0492 • 10" 5 


2.23 



Table 5: Difference between the gap-tooth scheme for (|30|) and a finite difference 
scheme for l|31|) for order k = 2 and Ax = 0.05 at time t = 2 • 10 -2 . 

h = 0.005 and 0.01. The results are shown in table 01 ^From this table we 
see that a smaller time-step indeed gives a smaller difference with respect 
to the corresponding finite difference scheme. However, we do not observe 
the ratio 2. We can show that this is due to the interference of the 0(-^)- 
term. Indeed, if we replace the microscopic problem inside each box with 
the homogenized problem with Neumann boundary conditions, the 0{-^)- 
term vanishes. The result is shown in table El We see that the decrease of 
convergence speed indeed disappears. The observed ratios are slightly larger 
than 2, due to the interference with the 0(h 2 ) term due to averaging, which 
is opposite in sign. 

5.3 Rough non-periodic (random) diffusion 

With this example, we illustrate that the scheme can also be used to simulate 
systems with rough coefficients, which are correlated on a length scale e. We 
use the same example as Abdulle and E pQ, who constructed it to demon- 
strate the behaviour of the heterogeneous multi-scale method. We first take 
a uniformly distributed random signal s(x) in [0.1, 1.1], with x E [0, 1]. We 
discretize the interval [0, 1] in N equidistant points Xi, and define a correla- 
tion kernel g e (x), such that 

g e (0) = -, e (x) = Oifxg(-^), [ 2 g e (x)6x = l. 
e 2 2 J_* 
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Figure 7: Left: Random correlated diffusion coefficient a e {x). Right: Reference 
solution (with N = 20001, At = 1 ■ 10~ 5 with an implicit Euler time-step) and 
gap-tooth solution (with Ax = O.f , At = 2 • 10 -3 ), at time t = 2 ■ 10~ 2 . 



with a e (x) constructed as above. Note that this is a deterministic problem 
once the diffusion coefficient is obtained. We consider only one realization 
of the diffusion coefficient field. Here, we can approximate the effective 
behaviour by averaging in space; in |20j . one had to average over many 
(shifted) initial conditions. 

We compute the solution with N = 20001 and e = 1 • 10~ 3 , and we com- 
pare this to a gap-tooth solution with h = 5 • 10~ 3 . As an initial condition, 
we take u e (x,0) = 1 — 4(x — 0.5) 2 , with homogeneous Dirichlet boundary 
conditions. Figure shows the diffusion coefficient and the reference and 
gap-tooth solutions. We compute the solution for this problem with an in- 
creasing number of boxes, where v = 0.2 is fixed and At = vAx 2 . The 
error is shown in table We see that the error decreases when more boxes 
are inserted. Note that, due to the roughness of the signal, it is difficult to 
draw conclusions on the convergence rate, and to determine good parameter 
values for the gap-tooth scheme. E.g. the length over which the gradient is 
averaged at the end points of each box is no longer uniquely defined, since 
the small scale e is only a correlation length. 



Here, we choose g e (x) = ^(1 — sin(^p)), x € (— f , §). We then define the 
rough diffusion coefficient in the discretization points as 



a e ( Xi ) = * s)(xi) = J\ g e ( Xi - £) S (£)d£. 



We then consider the diffusion equation 
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error 


Ax = 2 ■ 10 -1 
Ax = 1 • 10" 1 
Ax = 5 • 10~ 2 
Ax = 2.5 • 10" 2 


1.0283 • 10~ 2 
8.1064 • 10~ 3 
7.4622 • 10" 3 
1.4028 • 10" 3 



Table 6: Difference between the gap-tooth scheme with h = 5-10 3 and the reference 
solution at time t = 2 • 10~ 2 for increasing level of discretization. 

6 Avoiding the algebraic constraint 

We recall that the gap-tooth scheme, as presented above, performs the sim- 
ulations inside each box using an algebraic constraint, ensuring that the 
initial macroscopic gradient is preserved at the boundary of each box over 
the time-step At. If our goal is to accelerate time-integration using an ex- 
isting microscopic code, this constraint may require us to alter this code, so 
as to impose this macroscopically-inspired constraint. This may be imprac- 
tical (e.g. if the macroscopic gradient has to be estimated), undesirable (e.g. 
if the development of the code is expensive and time-consuming) or even 
impossible (e.g. if the microscopic code is a legacy code). 

Generally, a given microscopic code allows us to run with a set of pre- 
defined boundary conditions. It is highly non-trivial to impose macroscop- 
ically inspired boundary conditions on such microscopic codes, see e.g. |16| 
for a control-based strategy. This can be circumvented by introducing buffer 
regions at the boundary of each small box, which shield the dynamics within 
the computational domain of interest from boundary effects during a short 
time interval. One then uses the microscopic code with its built-in boundary 
conditions |2~T] . 

6.1 The gap-tooth scheme with buffers 

We note that, for a correct simulation, the only crucial issue is that the 
detailed system in each box should evolve as if it were embedded in a larger 
domain. This can be effectively accomplished by introducing a larger box 
of size H >> h around each macroscopic mesh point, but still only use 
(for macro-purposes) the evolution over the smaller, "inner" box. This is 
illustrated in figure |H1 Lifting and evolution (using arbitrary outer boundary 
conditions) are performed in the larger box; yet the restriction is done by 
taking the average of the solution over the inner, small box. The goal of 
the additional computational domains, the buffers, is to buffer the solution 
inside the small box from outer boundary effects. This can be accomplished 
over short enough times, provided the buffers are large enough; analyzing 
the method is tantamount to making these statements quantitative. 
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Figure 8: A schematic representation of the gap-tooth scheme with buffer boxes. 
The simulation is done in the box of size H, whereas for restriction, only information 
inside the inner box of size h is used. 

The idea of using a buffer region was also used in the multi-scale fi- 
nite element method (oversampling) of Hou ^U] to eliminate the boundary 
layer effect; also Hadjiconstantinou makes use of overlap regions to couple 
a particle simulator with a continuum code UJ. If the microscopic code al- 
lows a choice of different types of "outer" microscopic boundary conditions, 
selecting the size of the buffer may also depend on this choice. 

6.2 Damping factors 

Here, we show that we can study the gap-tooth scheme (with buffers) 
through its numerically obtained damping factors, i.e. by estimating the 
eigenvalues of its linearization. Integration with nearby coarse initial condi- 
tions is used to estimate matrix-vector products of the linearization of the 
coarse time-At map with known perturbation vectors; these are integrated 
in matrix-free iterative methods such as Arnoldi eigensolvers. For the dif- 
fusion homogenization problem (|30j) . we show that the eigenvalues of the 
gap-tooth scheme are approximately the same as those of the corresponding 
finite difference scheme for (|31[). When we impose Dirichlet boundary con- 
ditions at the boundary of the buffers, we show that the scheme converges 
to the standard gap-tooth scheme for increasing buffer size. 

Convergence results are typically established by proving consistency and 
stability. If one can prove that the error in each time step can be made 
arbitrarily small by refining the spatial and temporal mesh size, and that 
an error made at time t n does not get amplified in future time-steps, one 
has proved convergence. This requires the solution operator to be stable 
as well. In the absence of explicit formulas, one can examine the damping 
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factors of the time-stepper. If, for decreasing mesh sizes, all (finitely many) 
eigenvalues and eigenfunctions of the time-stepper converge to the dominant 
eigenvalues and eigenfunctions of the time evolution operator, one expects 
the solution of the scheme to converge to the true solution of the evolution 
problem. 

Consider equation (j31j) with Dirichlet boundary conditions u(0,t) = 
and u(l,t), and denote its solution at time t by the time evolution operator 

u(x,t) = s(u (x);t), (34) 

We know that 

s(sm(rrnrx); t) = e - ( m?r ) *sin(m7rx), m 6 N. 

Therefore, if we consider the time evolution operator over a fixed time i, 
s(-,t), then this operator has eigenfunctions sin(m-7rx), with resp. eigenvalues 

A m = e-t™) 2 *. (35) 

A good (finite difference) scheme approximates well all eigenvalues whose 
eigenfunctions can be represented on the given mesh. We choose t as a 
multiple of At for convenience. 

Since the operator defined in (|34[) is linear, the numerical time integration 
is equivalent to a matrix-vector product. Therefore, we can compute the 
eigenvalues using matrix-free linear algebra techniques, even for the gap- 
tooth scheme, for which it might not even be possible to obtain a closed 
expression for the matrix. Note that this idea is general; here we use it as 
a tool to study the effect of the buffer size on convergence of the gap-tooth 
scheme. However, although this analysis gives us an indication about the 
quality of the scheme, it is by no means a proof of convergence. 

6.3 Numerical results 

We illustrate this idea by computing the eigenvalues of the gap-tooth scheme 
of order k = 2, applied to (|3U|) . In this case, we know from section 0] that 
these eigenvalues should approximate the eigenvalues of a finite difference 
scheme on the same mesh. As method parameters, we choose Ax = 0.05, 
h = 5 TO" 3 , At = 2.5 T0~ 4 for a time horizon i = 4T0~ 3 , which corresponds 
to 16 gap-tooth steps. Inside each box, we use a finite difference scheme of 
order 2 with Sx = 1 ■ 10 -5 and an implicit Euler time-step of 5 • 10~ 7 . 
We compare these eigenvalues to those of the finite difference scheme with 
Ax = 0.05 and At = 2.5 • 10~ 4 , and with the dominant eigenvalues of the 
reference solution (a finite difference approximation with Ax = 1 • 10 -5 , 
At = 5 • 10 -7 and implicit Euler time-stepping). The result is shown in 
figurelHl We now introduce a buffer region of size H, and we impose Dirichlet 
boundary conditions at the outer boundary of the buffer region. Lifting is 
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Figure 9: Comparison between the damping factors (left) and the eigenfunction 
corresponding to eigenvalue A 3 (right) of the exact solution (full line), the finite 
difference approximation (dashed) and the gap-tooth scheme (dotted). The eigen- 
function of the gap-tooth scheme is indistinguishable of the finite difference eigen- 
function. 

done in identically the same way as for the gap-tooth scheme without buffers; 
we only use © as the initial condition in the larger box [xi — y, X{ + y]. We 
compare the eigenvalues again with the equivalent finite difference scheme 
and the exact solution, for increasing sizes of the buffer box H. Figure El] 
shows that, as H increases, the eigenvalues of the scheme converge to those 
of the original gap-tooth scheme. We see that, in this case, we would need 
a buffer of size H = 4 • 10 -2 , i.e. 80% of the original domain, for a good 
approximation of the damping factors. It is possible to decrease the buffer 
size by decreasing At, which results in more re-initializations. 

7 Conclusions 

We described the gap-tooth scheme for the numerical simulation of multi- 
scale problems. This scheme simulates the macroscopic behaviour over a 
macroscopic domain when only a microscopic model is explicitly available. 
We analyzed the convergence of this scheme for a parabolic homogenization 
problem with non-linear reaction. 

We showed that our method approximates a finite difference scheme of 
arbitrary (even) order for the homogenized equation when we appropriately 
constrain the microscopic problem in the boxes, and illustrated this theo- 
retical result with numerical tests on several model problems. Our analysis 
revealed that the presence of microscopic scales, combined with the require- 
ment that the macroscopic gradient does not change over one gap-tooth 
time-step At, introduces an error term that grows with decreasing At, which 
is not optimal. 

We also demonstrated that it is possible to obtain a convergent scheme 
without constraining the microscopic code, by introducing buffers that shield 
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Figure 10: Comparison between the damping factors (left) and the eigenfunction 
corresponding to the eigenvalue A 3 (right) of the exact solution (full line), the finite 
difference scheme (dashed) and the gap-tooth scheme with buffers (dash-dotted 
lines) for increasing buffer sizes H = 2 • 10~ 2 , 3 • 10~ 2 . . . , 1 • 10 _1 . 

over relatively short time intervals the dynamics inside each box from bound- 
ary effects. It is possible, even without analytic formulas, to study the 
properties of the gap-tooth scheme and generalizations through the damp- 
ing factors of the resulting coarse time-Ai map. In a forthcoming paper, 
we will use these damping factors to study the the trade-off between the 
effort required to impose a particular type of boundary conditions and the 
efficiency gain due to smaller buffer sizes and/or longer possible time-steps 
before reinitialization. 

The time-stepper as constructed in this paper will allow us to perform 
simulations of the effective behaviour of a microscopic system over macro- 
scopic space and macroscopic time (when combined with projective integra- 
tion), or to perform tasks as bifurcation analysis or coarse control (when 
coupled to time-stepper based bifurcation codes). 
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